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ABSTRACT 

We study the evolution of phase-space density during the hierarchical structure formation 
of ACDM halos. We compute both a spherically-averaged surrogate for phase-space density 
(Q = p/'^^) ™d the coarse-grained distribution function /(x, v) for dark matter particles that 
lie within ^ 2 virial radii of four Milky-Way-sized dark matter halos. The estimated /(x, v) 
spans over four decades at any radius. Dark matter particles that end up within two virial 
radii of a Milky-Way-sized DM halo at z = have an approximately Gaussian distribution in 
log(/) at early redshifts, but the distribution becomes increasingly skewed at lower redshifts. 
The value /peak corresponding to the peak of the Gaussian decreases as the evolution pro- 
gresses and is well described by /pcak(-z) oc (1 + z)'^-'^^^-^. The decrease is due both to the 
dynamical mixing as the matter accreted by halos is virialized, and due to the overall decrease 
of space-density of the unprocessed material due to expansion of the Universe. The highest 
values of /, (responsible for the skewness of the profile) are found at the centers of dark matter 
halos and subhalos, where / can be an order of magnitude higher than in the center of the main 
halo. We confirm that Q (r ) can be described by a power-law with the slope of/3 = — 1.8±0.1 
over 2.5 orders of magnitude in radius and over a wide range of redshifts. This Q{r) profile 

likely reflects the distribution of entropy {K = cr^/p^^ oc r^-^), which dark matter acquires 
as it is accreted onto a growing halo. The estimated /(x, v), on the other hand, exhibits a more 
complicated behavior Although the median coarse-grained phase-space density profile F{r) 
can be approximated by a power-law, oc ^--i 6±o.i5^ jj^g inner regions of halos (< 0.6 Tvir), 
at larger radii the profile flattens significantly. This is because phase-space density averaged 
on small scales is sensitive to the high-/ material associated with surviving subhalos, as well 
as relatively unmixed streams resulting from disrupted subhalos, which constitute a sizeable 
fraction of matter at large radii. 

1 INTRODUCTION tributions of cosm ological DM halos are described by si mple 3- 

parameter profiles jNavarro et al.lll996l : lMerritt et alJlOOd) . 

Cosmological TV-body simulations of the formation of structure in The study of the evolution of the phase-space distribution 

the Universe allow us to reconstruct how gravitationally bound ob- function (DF) of a collisionless dynamical system (composed of 

jects like galaxies and clusters form and evolve. These simulations stars or DM or both) is fundamental to understanding its dynami- 

have shown that, despite the seemingly complex hierarchical for- cal evolution, since a collisionless system is completely described 

mation history of dark matter (DM) halos, the matter density dis- by its phase-space density distribution (otherwise called the fine- 
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grained DF f (x, v)). For an isolated collisionless system, f (x, v, t) 
is fully described by the Collisionless Boltzmann Equation (CBE) 
(sometimes called the Vlasov equation), which is a continuity equa- 
tion in phase-space. A consequence of the CBE is that the fine- 
grained DF f is always conserved. In addition, V(f )cif , the volume 
of phase-space occupied by phase-space elements whose density 
lies between (f, f + df), is also conserved. However, the conser- 
vation of the fine-grained f (x, v) and V(f)df is not a very useful 
property for understanding the evolution of real systems, since in 
practice it is only possible to compute an average of f over a finite 
volume of phase-space. This average is referred to as the coarse- 
grained DF /(x, v), and the associated volume density is V{f). 
The evolu tion of the coarse-grained DF is governed hy Mixing 
Theorems teinnev & TremainejigSTl i lTremaine et al.ll986l : lMathutl 
Il988l) . which state that if the coarse-grained DF is bounded, then 
processes that operate during the relaxation of collisionless systems 
(e.g. phase mixing, chaotic mixing, and the mixing of energy and 
angular momentum that accompanies violent relaxation) result in a 
decrease in all convex functions of the coarse-grained phase-space 
density. In particular, the coarse-grained phase-space density in a 
region of phase-space is expected to decrease with time, because at 
any point in phase-space , both low and high phase-space density re- 
gions get highly mixed ( Binnev & Tremairigll987ljTremaine et al.l 
Il986h . In other words, mixing reduces / such that, at any time, 
/ < /max, the maximum phase-space density at the start. In ad- 
dition, there are simple relationships that exist between the initial 
V(f) and the coar se-grained Vi f) at a later time: V{f ) is always 
greater than V(f ) ( lMathuij[l98a) . 

Recently, 'De hnerj ( l2005h stated and proved a new Mixing 
Theorem that states that the excess-mass function D{f), which 
is the mass of material with coarse-grained phase-space density 
greater than a value /, always decreases due to mixing, for all 
values of /. The excess-mass function is additive so that the ex- 
cess mass of a combination of non-overlapping systems is the sum 
of their individual excess masses. He showed using this theorem 
that steeper cusps are less mixed than shallower ones, independent 
of the details of the DF or density profile, and this implies that a 
merger remnant cannot have a cusp that is either steeper or shal- 
lower than the steepest of its progenitors. This theorem has pow- 
erful implications for the evolution of halos with cuspy power-law 
central phase-space density profiles via mergers, once cuspy pro- 
files have been formed. 

iTavlor & Navarrol ( I2OO1I) have introduced a surrogate quantity, 
Q{r) = p{r)/a{r)^ , where p(r) is the configuration-space density 
averaged in spherical shells, and a{r) is the velocity dispersion of 
dark matter particles averaged in a spherical shell centered at the 
radius r. Note that, although Q{r) has dimensions of phase-space 
density, it is definitely not true coarse-grained phase-space density, 
because it is constructed out of separately computed configuration- 
space density and velocity dispersion in arbitrarily chosen spher- 
ical shells. Consequently, it is not expe cted to obey any Mixing 
Theorems ( iDehnen & McLaughlinlboOSi) and is generally difficult 



to interpret in terms of phas e-space density. Neverthel ess, it is an 
intriguing quantity, because iTavlor & Navarrol i200lh show that 
Q(r) is well approximated by a single power-law, Q oc r"'' with 
(3 ~ 1.874 over more than 2.5 orders of magnitude in radius for 
CDM halos univ ersally, regardless of mass and backgro und cos- 
mology (see also lRasia et al.ll2004l ; lAscasibar et al.ll2004) . Recent 
work, however, indicates that Q{r) profiles are somewhat sensitive 
to the amount of subst ructure and the slope o f the power spectrum 
dWang & White 2008i ; lKnollmann et ai]|2008l) . 

While power-law Q{r) profiles are by no me ans universal 
for se lf-gravitating systems in dynamical equilibrium dBames et al.l 
12006*), power-law profiles with exactly the same slope were first 
shown t o arise in simple se lf-similar spherical gravitational infall 
models jBertschingeij 1 1 985h pointin g to a possible univ ersality in 



the mechanisms th at produces them. [ Austin et extended 
the early work of iBertschingeJ ( 1 19851) using semi-analytical ex- 
tended secondary infall models to follow the evolution of collision- 
less spherical shells of matter that are initially set to be out of dy- 
namical equilibrium and are allowed to move only radially. They 
concluded that the power-law behavior of the final phase-space 
density profile is a robust feature of virialized halos, which have 
reached equilibrium via violent relaxation. They used a constrained 
Jeans equation analysis to show that this equation has different 
types of solutions and that it admits a unique periodic solution 
which gives a power-law density profile with slope /3 — 1.9444 
(comparable to the num erical value of P — 1.87 obtained by 
ITavlor & Navarrol ( 1200 ih ). 



A recent study of lHoffman et al.l J2007h followed the evolution 
of Q{r) profiles in cosmological DM halos in constrained simu- 
lations designed to control the merging history of a given halo. 
These authors showed that, during relatively quiescent phases of 
halo evolution, the density profile closely follows that of a NEW 
halo, and Q{r) is always well represented by a power-law r"'' 
with P — 1.9 ± 0.1. They showed that Q(r) deviates from power- 
law most strongly during major mergers but recovers the power-law 
form thereafter (this is consistent with our findin gs using a series 
of controlled merger simulations, I Vass et al]|2008i) . More recently, 
jWoitak et al.l 1200^ demonstrated (using a varient of the FiEs- 
tAS code) that the DF of ACDM halos of mass lO" - 10^^ Mq can 
be separated into energy and angular momentum components and 
proposed a phenomenological model for spherical potentials. They 
showed that their model DF was a good match to the N-body DEs 
and reproduced the power-l aw behavior of Q(r). Another stud y 
of relevance is the work of IPeirani & de Freitas Pachecd ( 120071) . 
which defined a global value for the phase-space density, Q, which 
they found decreases rapidly with time. Despite the insights ob- 
tained in these studies, the origin for such universality of Q{r) is 
not yet understood, and the question of how this quantity relates to 
the true coarse-grained phase-space density has not yet been thor- 
oughly investigated, which motivates our study in which we focus 
on the evolution of both Q{r) and phase-space density evolution 
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using a set of self-consistent cosmological simulations of halo for- 
mation. 

In the last four years, three independent numerical codes 
to compute the true coarse-grained phase-sp ace density from 

., dA 

lAscasibar & BinnevI l2005l : Isharma & Steinmetj |2006|) . The first 



A^-body simulations hav e been developed 



Arad et al. 



y trom 
f |2004 



study of coarse-grained / showed that it has its highest values at the 
centers of DM halos and subhalos. They also showed that the vol- 
ume density of phase-space V{f) within individual cosmological 
dark matter halos at 2 = has a power-law p rofile over nearly four 
orders of magnitude in / ( lAradetalj|2004h . The main criticism 
of this code (which is based on the Delunay tesselation of phase- 
space) is that it is not metric free. In this paper, we compare results 
obtain ed with two o ther codes "FiEstAS " l|Ascasibar & BinnevI 
l2005h and "EnBiD " jSharma & Steinmetj|2006l) , which are both 
fast and metric fre^E to study the evolution of the course-grained 
phase-space density. 

In an accompanying paper jVass et al.l2008h . we present anal- 
ysis of phase-space and Q{r) evolution during major mergers be- 
tween CDM-like spherical halos. We show that major mergers 
(those which are the most violent and therefore likely to result in 
the greatest amount of mixing in phase-space) preserve spherically- 
averaged phase-space density profiles (both Q{r) and F{r)) out to 
the virial radius r^h of the final remnant. Mergers between identical 
but non-NFW halos, i.e. halos with central cusps significantly shal- 
lower than NEW (7 = 0.2) or steeper than NEW (7 = 1.7), also 
retain their original phase-space density profiles. Mergers between 
a halo with a shallow cusp and a steep cusp produce halos with in- 
ner phase-space density profiles that are indistinguishable from that 
of the steepest cusp in the pair of p rogenitors. These results were 
anticipated by the work of lDehned (2005) and are a consequence 
of the Mixing Theorem for cusps, which requires that the excess 
mass function of collisi onless systems always decreases. Thus we 
confirm the prediction toehnerj ( |2005|) ) that the phase-space den- 
sity profiles of DM halos are extremely robust, and in particular, 
the steepest central cusp always survives. 

To understand th e robustness of profiles outsi de the center 
cusp Vass et all ^20081) , we drew on recent work by IValluri et al.l 
1200/) that shows that matter in equally-spaced radial shells is re- 
distributed during the merger is such a way that only about 20% of 
the matter in the central cusp is ejected during the merger to radii 
extending out to about two scale radii. Eor all other radii, roughly 
15% by mass in each radial shell is redistributed uniformly with 



^ An ev en more sophisticated nu merical method has recently been pre- 
sented bv'Vogelsberg er et al.l j2008l) for calculating the fine-grained phase- 
space structure of dai'k matter distributions derived from cosmological sim- 
ulations. This code has the potential to identify fine-scale structure such as 
caustics in phase-space and the phase-space structure of tidal streams in the 
Milky- Way halo. Its ability to estimate the fine-grained f(x,v) makes it 
useful for understanding non-equilibrium systems and for making accurate 
predictions for direct DM searches. 



radius. Eurther, nearly 40% of the mass of the pr ogenitor halos lies 
beyond the formal virial radius of the remnant JKazantzidis et al.l 
l200dlVaUurietalJl2007h . and this matter originates roughly uni- 
formly from each radial interval starting from about three scale 
radii from the center . While major me rgers, such as those stud- 
ied bv lDehnerJ j2005l) : IVass et al.rj2008h . are instructive, since they 
represent the most extreme form of mixing, they are not the major 
mode of mass accretion in the Universe. In addition, since experi- 
ments with major mergers demonstrate that pre-existing power-law 
profiles are preserved, they shed little light on the formation of the 
profiles. Our primary goal is to gain a better understanding of the 
origin of the power-law phase-space density profiles seen in cos- 
mological DM simulations. 

In this paper we investigate the evolution of the coarse-grained 
phase- space density in the formation and evolution of four Milky- 
Way-sized halos in a ACDM cosmology with cosmological param- 
eters: {Qrn,^A,h,as) ~ (0.3, 0.7, 0.7, 0.9). In § |2] we summa- 
rize the cosmological N-body simulations analyzed in this study as 
well as the numerical methods used to obtain coarse-grained phase- 
space densities. (A detailed comparison of the two codes used (Ei- 
EstASand EnBiD) is presented in the Appendix, where we give 
reasons for our choice of EnBiD , with a n = 10 smoothing ker- 
nel). In §[3] we describe the properties of the spherically-averaged 
quantity Q(r) as well as the the coarse-grained DE for one Milky- 
Way-sized DM halo al z = 0, as well as the evolution of the phase- 
space DE of this halo from 2: = 9 to the present. In §|4]we compare 
these results for three additional Milky- Way- sized halos. In §|5]we 
present an interpretation of the power-law profiles seen in Q{r) 
and discuss the implications of observed coarse grained distribu- 
tion function in a cosmological context. §[6] summarizes the results 
of this paper and concludes. 



2 NUMERICAL METHODS 

The simulations an alyzed in this paper ar e des cribed in greater de- 
tail in the works of iKravtsov et al.l j2004) and lGnedin & Kravtsovl 
(f200§). The simulations were carrie d out using the Adapti ve Re- 
finement Tree iV-body code (ART, IKravtsov etafl ( Il997l) ). The 
simulation starts with a uniform 256'^ grid covering the entire com- 
putational box. This grid defines the lowest (zeroth) level of resolu- 
tion. Higher force resolution is achieved in the regions correspond- 
ing to collapsing structures by recursive refining of all such regions 
by using an adaptive refinement algorithm. Each cell can be refined 
or de-refined individually. The cells are refined if the particle mass 
contained within them exceeds a certain specified value. The grid is 
thus refined to follow the collapsing objects in a quasi-Lagrangian 
fashion. 

Three of the galactic DM halos were simulated in a comoving 
box of 25 h^^ Mpc (hereafter L25); they were selected to reside 
in a well-defined filament at 2 = 0. Two halos are neighbors, lo- 
cated 425 h^^ kpc from each other. The third halo is isolated and 
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is located 2 Mpc away from the pair. Hereafter, we refer to the iso- 
lated halo as Gl and the halos in the pair as G2 and G3. The virial 
m asses and virial radii for the halos studied are given in Table 1 
in lKravtsov et alj j2004) . The virial radius (and the corresponding 
virial mass) was chosen as the radius encompassing a mean density 
of ~ 200 times the mean density of the Universe. The masses of 
the DM halos are well within the range of poss ibilities allowed b y 
models for the halo of the Milky Way galaxy (Klvpin et al. 2002). 
The simulations followed a Lagrangian region corresponding to the 
sphere of radius equal to two virial radii around each halo at z = 0. 
This region was re-sampled with the highest resolution particles of 
mass TTip — 1.2 x 10''/i~^Mq, corresponding to 1024'' particles 
in the box, at the initial redshift of the simulation {zi — 50). The 
maximum of ten refinement levels was reached in the simulations 
corresponding to the peak formal spatial resolution of 150 comov- 
ing parsec. Each host halo is resolved with ~ 10® particles within 
its virial radius at 2; = 0. 

The fourth halo was simulated in a comoving box of 20 h^^ 
Mpc box (hereafter L20) and it was used to follow the Lagrangian 
region of approximately five virial radii around the Milky Way- 
sized halo with high resolution. In the high-resolution region the 
mass of the dark matter particles is mp = 6.1 x 10^ /i^^ Mq, 
corresponding to effective 1024^ particles in the box, at the initial 
redshift of the simulation (21 — 70). As in the other simulation, 
this run starts with a uniform 256^ grid covering the entire com- 
putational box. Higher force resolution is achieved in the regions 
corresponding to collapsing structures by recursive refining of all 
such regions using an adaptive refinement algorithm. Only regions 
containing highest resolution particles were adaptively refined. The 
maximum of nine refinement levels was reached in the simulation 
corresponding to the peak formal spatial resolution of 150h~^ co- 
moving parsec. The Milky Way-sized host halo has the virial mass 
of 1.4 X 10^^ h^^ Mq (or 2.3 million particles within the virial 
radius) and virial radius of 230/i~^kpc at 2 = 0. 

The peak spatial resolution of the simulations determines the 
minimum radius to which we can trust the density and velocity pro- 
files. In the following analysis, we only consider profiles at radii at 
least eight times larger than the peak resolution of the simulations 
(i.e., « 0.5 — Q.6h~^ comoving kpc or r ~ 0.004 rvir(z=o))- 



2.1 Phase-space density estimators 

In this paper, we will focus on the time evolution of the spher- 
ically averaged phase-space density Q{r) as well as the coarse- 
grained phase-space DF /(x, v) during the formation of the Milky- 
Way sized DM halos in ACDM cosmological simulations described 
above. 

In the last four years three independent numerical codes 
to compute the coarse-grained phase-space density from A^- 
body simulations have b een developed 1 Arad et al] |2004| : 
lAscasibar & BinnevI |2005| ; Isharma & Steinmetj I2OO6I) . These 



techniques differ in the scheme used to tesselate 6-dimensional 
phase space as well as in the density estima tors they use. In our 
analysis, we use the FiEstAS algorithm of Ascasibar & Binne^ 



(2005) and the EnBid algorithm of Sharma & Steinmetz (200i 
due to their speed and the metric free nature. The FiEstAS 
algorithm is based on a repeated division of each dimension of 
phase-space into two regions that contain rough ly equal numbers 
of particles. EnBiD (Shar ma & Steinmetz! [2006) closely follows 
the method used by FiEstAS , but improves upon it to allow for 
more accurate computation of / at high phase-space densities. 
In EnBiD the space is tesselated along the dimension having 
minimum entropy (and maximum information) each step of the 
computation. This scheme optimizes the number of divisions to be 
made in a particular dimension and extracts maximum information 
from the data by using a minimum entropy criterion based on the 
Shannon entropy that achieves much greater accuracy a t measuririg 
the phase-space density when it is high. ISharma & Steinmetg 

(2006) showed that EnBiD with a kernel which includes 10 nearest 
neighbors (n = 10) about each point was able to recover analytic 
phase-space density profiles to nearly 3-4 decades higher values of 
/ than FiEstAS . 

For our halos at 2 = 0, the results we obtained with the dif- 
ferent codes were complet ely consistent with those obtained by 
[sii arma & Steinmetz! ( l2006h . The deviations between the estimates 
obtained with the different codes and various parameters were sig- 
nificantly higher at higher redshift - where comparisons with an- 
alytic estimates are not available. We refer the reader to the Ap- 
pendix for a detailed comparison between estimated values of / 
using the two codes at various redshifts. In the rest of this paper we 
present results obtained with EnBiD (n = 1 kernel), the parame- 
ters prefered bv lSharma & Steinmetz! ( 12006!) . 

Figure \T\ Top shows the spherically averaged quantity Q — 
p/u"^ (dot-dashed curve); the mean of /(x,v) in spherical bins 
(dashed curve) and the median DF in spherical bins (dotted curve). 
As we will show in Figure |4l the large fluctuations in the mean 
value of / (dashed curve) beyond O.lr-vir are due to the presence 
of substructure, which can have extremely high central values of 
/. The median value of / (hereafter represented by F) is much 
smoother, being less sensitive to the large range (nearly eight or- 
ders of magnitude) in / at each radius. In what follows, we will 
use the median F, computed in concentric radial bins centered on 
the most bound particle in the main halo, since it is less sensitive to 
substructure. Q{r) is well fitted by a power-law: Q oc ^-i-*4±o.oi2 
over the radial range [0.004-0.6] r/r^ii (thin red line). A power- 



law is also fitted to F{r) oc r 



'1.59±0.054 



(blue line) over the same 



radial range. The bottom panel shows the residuals of the fits F 
(log(F/Ffit)) and Q (log(Q/Qfit))- This plot shows that while 
Q(r) is an extremely good power-law, the median F{r) is only ap- 
proximately power-law over the same radial range. A s imilar result 
was obtained in a much higher resolution simulation by Sta del et al.l 
[2008.) . a study which appeared as we were preparing this paper for 
publication. 
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Figure 1. Q{r) and F{r) in 100 spherical radial bins about the most- 
bound-pai1icle in halo Gl from the L25 simulation at z=0. Top panel: Q 
(dot-dashed curve). Power-law fit Q oc r~^'^'^^'^''^^^ is given by solid red 
line. F(r) from EnBiD (with n = 10 smoothing kernel (dotted curve)). 
Mean values of f estimated by EnBiD (n = 10)are given by long-dashed 
curve is noiser due to subhalos. A power-law fit to F{r) oc r-l-59±0.054 
is given by solid blue line. Residuals to the two power-law fits are shown in 
the bottom panel. 



EVOLUTION OF PHASE-SPACE DENSITY WITH 
REDSHIFT 



studied the phase-space density profiles of a 
dark matter halo by tracking Q{r) for material within the formal 
virial radius at each redshift. They found that the virialized mate- 
rial within this radius has an approximately power-law from with 
a constant power-law index of —1.9 ± 0.05 at all redshifts from 
= 5 to the present. 

In this paper we follow a slightly different approach since our 
objective is to understand the evolution of the true coarse-grained 
phase-space density distribution. We track the evolution of phase- 
space density, by tracing backwards in time, all the material that lies 
inside the virial radius at 2 = 0. This will allow us to understand 
how the initially high phase-space density material that lies out- 
side virialized systems at high redshift falls in and undergoes mix- 
ing and how the resultant mixing preserves the phase-space density 
profiles as a function of redshift. 

We identified Milky- Way sized DM halos at 2 = in our cos- 
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Figure 2. as a function of co-moving radius r for all material that lies 
within iTvir of the center of the Gl halo at 2 = 0. Each curve coiTesponds 
to a different redshift as indicated in the legends. 



mological simulations and identify all the particles that lie inside 
twice the virial radius at 2 = 0. We then tracked these particles 
back to 2 — 9. Our simulations do proceed back in time to even 
higher redshifts, however we do not analyze them here because the 
mass resolution at higher redshifts does not allow us to resolve the 
evolution of most of the objects beyond this epoch. In total there 
were 1.6x 10® particles within 2 virial radii of the Gl halo, which 
is the halo that is most isolated at 2 = 0. Two other halos (G2 and 
G3) from the L25 simulation, are closer to each other at 2 = 
(1.5 virial radii apart), and thus, we only track back particles within 
one virial radius of the centers of these two halos. The particles in 
the high resolution simulation L20 were selected and tracked in the 
same way as for the Gl halo, and at 2 = there were 3 x lO*" par- 
ticles inside twice the virial radius. The position and velocity data 
for all the particles identified as belonging to a given halo at 2 = 
were analysed in physical coordinates (distances are always in kpc 
and velocities in kms^^). In this section, we will present the results 
obtained for the Gl halo in the L25 simulation. A comparison with 
results for the other three halos is made in §|4] 

Figure [2] plots the median value (F) of the phase space den- 
sity / as a function of distance to center of the most massive pro- 
genitor at several different redshifts. At the highest redshift plotted 
(2 — 9) the overall value of F is higher at large radii than it is at 
the center of the halo and only varies by a factor of few over the 
entire range of radii (i^ ~ 10 — 3OAf0kpc^'^(kms~^)~^ ). This 
is because within the inner, virialized regions the coarse-grained 
phase space density is lowered due to mixing, while it is still high 
in the outer regions where the large fraction of matter is not yet 
mixed or is located in small subhalos, which have not mixed to the 
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degree of the main host, and the tidal streams they produce when 
disrupted. As the evolution progresses, there is a decrease in the 
central value of F till « = 3. The median F at the center of the 
halo drops from F « lOM0kpc~"^(kms~^)~"^ at 2: = 9 down 
to F f» 2Mokpc"^(kms"^)"^ at z = 0. Beyond z = 3, the 
central value remains constant while there is a steady decrease in 
F at larger radii leading to the development of an approximately 
power-law profile at 2: = 0. 

Although it is instructive to plot the mean and/or median val- 
ues of F at each radius, much more information is contained in 
the full coarse-grained phase-space density /. The six-dimensional 
function /(x, v) can be most easily visualized us ing the volum e 
DP y(/) which also obeys a Mixing Theorem jMathuJ IT988h . 
lArad et d] fe004) showed that at z = 0, V{f) for cosmological 
dark matter halos follows a power-law profile V{f) oc /~° with 
power-law index a — 2.4 over 4 orders of magnitude in /. They ar- 
gued that since DM halos are almost spherical, if their phase-space 
DFs are approximately isotropic their DFs could be written as func- 
tions of energy alone: / — f{E). In this case they showed that if 
Q oc r"*^, then V{f) would also be described by a power-law, and 
the index a was related to the index /3 through a simple equation. 

In Figure [3] we plot the evolution of V{f) with redshift for 
all the material that lies within twice the virial radius of halo Gl 
at z = 0. In the bottom right hand panel, the solid line is a fit to 
V{f) for values of 10"" < / < 10^ ^ Over this range the best fit 
power-law profile is given by !/(/) cx; j -^-34±o.02 close 
to the power-law profile fit obtained bv lArad et alj ( l2004h . At all 
redshifts V{f) deviates from the power-law profile at both low and 
high end. It is likely that at the high-/ end the distribution deviates 
from the power-law due to the unresolved subhalos below the mass 
resolution limit of the simulation. 

The fine-grained version of the volume V(f) of phase-space 
associated with material of phase-space density f is conserved by 
the collisionless Boltzmann equation, however the coarse-grained 
version of this quantity (V{f)) always increases by the Mixing 
Theorem (Mathur I988|). We see that increase in V{f) occurs quite 
steadily from 2 = 9 onward at all values of /, but particularly for 
lower values of /. As the system evolves there is more volume as- 
sociated with material with low /. At any value of / the volume 
V{f) increases with decreasing redshift. Thus the behavior of the 
evolution of / with redshift is completely consistent with our ex- 
pectations from the Mixing Theorems. While the Mixing Theorems 
have previously been demonstrated for isolated systems, this is to 
our knowledge the first demonstration of their validity in a cosmo- 
logical context. 

It is particularly illuminating to plot the full coarse-grained 
phase space density of all particles as a function of radius from 
the center of mass of the main halo at each redshift (Figure |4l(. At 
any radius from the center there exist particles with phase-space 
densities spanning between 4-8 decades in /. The colored con- 
tours represent a constant particle number per unit area on the plane 
(log(/), log(r)). The yellow/orange contours represent the param- 



eter range with the largest number of particles while the blue/black 
contours represent regions with the smallest number of particles. 
The solid white curves on each plot show the median value (F) of 
/(x, v) in 100 logarithmically spaced bins in r. These curves are 
identical to the various curves in Figure|2]and trace the regions with 
the largest particle number density at each radius. Numerous spikes 
in / are seen at r > O.lrvir and correspond to DM sub-halos. 

As the evolution proceeds, there is an overall lowering of the 
median phase-space density F (white curves). There is also a con- 
tinuous decrease of the low-end envelope of / to lower values 
pointing to an increasingly large amount of matter with low phase- 
space density. However the upper envelope representing the highest 
phase space density regions at the centers of the subhalos remains 
relatively unchanged with redshift indicating that primordial val- 
ues of / are largely preserved at the centers of DM subhalos that 
lie outside the center of the main halo. 

It is illustrative to compare the evolution of F(r) and Q{r) 
with redshift. Figure [5] shows the two curves plotted as a function 
radius (in co-moving units) at each redshift for all the matter in 
halo Gl that lies within two virial radii at z = 0, as a function 
of physical radius from the center of the most massive progenitor 
of the final halo (in units of the virial radius at 2 = 0). In each 
panel a thin vertical line is drawn at the virial radius of the halo 
at that redshift. We note that as expected from the expansion of 
the Universe and increase in clustering, the virial radius of the halo 
increases with decreasing redshift. 

We find the best fit power-law Q oc r"^_f^^^;^]f_at2 = t o 
the profile within 0.6rvir. In agreement with lHoffman et al] ( l2007h . 
we find that the same power-law provides a reasonably good fit to 
the profiles at redshifts from z ~ 5 to 2 = 0. While Q always 
decreases with radius since it is a spatial average over increasingly 
large volumes of configuration space, the median phase-space den- 
sity F decreases monotonically with radius only within the virial 
radius of the main progenitor at that redshift. Outside the virial ra- 
dius at each redshift (i.e. to the right of the vertical dashed line), 
F{r) become significantly flatter, even increasing with increasing 
radius. At all 2 7^ this represents the median phase-space density 
of matter that will lie within ~ 1 virial radius at 2 = 0, but is either 
still unvirialized or lies within small halos which have not mixed 
to the degree of the more massive main progenitor. The nearly con- 
stant value of F beyond the virial radius at high 2 may be inter- 
preted as the median phase space density of DM in the Universe 
at that redshift. At lower redshifts, a non-negligible fraction of this 
matter lies within inner regions of small halos, which have under- 
gone relatively small amount of mixing. As more and more material 
undergoes substantial mixing as the halos evolve, the high phase- 
space density unmixed material in the subhalo centers, as well as 
in the relatively unmixed streams formed from disrupted subhalos, 
prevents the median at large radii from decreasing rapidly, despite 
the large increase in low material with low /-values within the viri- 
alized regions of the main halo. 

To better understand the evolution of / with redshift we plot 
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Figure 3. The volume DF V{f), for different redshifts in the evolution of the Gl halo in the L25 run. A power- law fit V{f) oc /-2.34±0.017^ jjje best fit 
over this range of /, is shown at 2: = (solid line). 



histograms of log(/) at each redshift in Figure [6] (solid line). The 
initial distribution of log(/) at 2; = 9 is close to Gaussian (indi- 
cating that / itself has a log-normal distribution). The peak of the 
Gaussian (the median of log(/)) moves to lower values as the halo 
evolves from z — 9 to z — due to both mixing in the main halo 
and in subhalos and due to cosmic expansion. 

By approximately z — 1, the skewness of the distribution is 
significant because more and more matter undergoes mixing dur- 
ing halo mergers and moves to lower values of /, and the skewness 



grows steadily until z — 0. The figure shows, however, that some 
fraction of the high-/ matter survives to z = 0. The thick dashed 
line shows the distribution of dark matter particles that have the 
values of / ^ lAfQkpc^^(kms~^)~^ (situated in the high / tail 
of the histogram) at z = 0. The thin blue dashed line plots the same 
distribution but with y-axis values multiplied by 10 to make the dis- 
tribution more visible. As we saw in Figure|4l most of these high-/ 
tail particles are in subhalos at z = 0, while some are also in the 
highest phase-space density particles in the central cusp of the halo 
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Figure 4. Contours of constant particle number density in the plane (log(/), log(r)) for the Gl halo, plotted at different redshifts. At each redshift the radial 
distribution of particles is given relative to the most bound halo pailicle in units of the virial radius of the halo z = 0. The yellow and black contours 
correspond to the regions with the highest and lowest particle number densities respectively. Contours are spaced at logarithmic intervals in particle number 
density relative to the maximum density contour. F, the median value of /(x, v) is represented by the solid white line. 



at 2 = 0. This high-/ sub-population has a Gaussian distribution at 
2 = 9 with a mean / = 1.73 ± O.6M0kpc"^(kms"^)'^ com- 
pared with the mean / = 1.57 ± O.54M0kpc"^(knis"^)"^ for 
all the particles (in the solid curve) at 2 = 9. A Students' T-test in- 
dicated with 99% confidence that both distributions are drawn from 



the same population at z — 9. This implies that the material that 
lies in the centers of DM subhalos at z = will have phase-space 
densities that are representative of the mean phase-space density of 
matter at z = 9. We expect that at even higher redshifts, both distri- 
butions in log(/) will become much more sharply peaked asymp- 
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Figure 5. Q = p/cr^ (dot-dashed line) and median F (dotted line) as a function of r lr^\^{z = 0) for different redshifts for the Gl halo. The dashed vertical 
line is situated at the virial radius of the virialized main halo at each redshift. 



toting to a narrow function at the start of the simulation, whose 
width represents the true range of / for CDM particles. 

This decrease in F of "unprocessed" material in Figure[5]and 
the decrease (see Fig.|6ll of peak of log(/) is a result of two inde- 
pendent factors. First, mixing (both due to merging between sub- 
halos and due to merging of subhalos with the main halo), causes 
a decrease in the overall coarse-grained phase-space density, as an- 
ticipated from the Mixing Theorem. Second, at most redshifts, a 
significant fraction of matter has not yet been virialized or is witin 
much smaller halos and lies outside the virial radius at that red- 



shift (to the right of the vertical line in Fig. [5). The density of 
this unprocessed material decreases due to the cosmic expansion as 
(.p{z) oc (1 + Thus the decrease in median phase-space den- 

sity of a halo with redshift is the consequence of both relaxation 
and mixing within virialized halos and subhalos, and the cosmic 
expansion of the Universe. 

^ Hubble expansion does not contribute to a change in the local velocity 
dispersion of dark matter particles that contributes to phase space density. 
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Figure 6. Histograms of the phase-space density / of all halo pailicles as a function of redshift (solid black curves). The thick blue dashed curves follow 
those pailicles which have the highest phase-space densities / ^ 1 at 2 = 0. (The thin blue dashed hne is the same as the thick-dashed line, except that the 
corresponding j;-axis values multiplied by 10 to enhance visibility.) 



4 COMPARSION OF FOUR MILKY- WAY SIZED HALOS 

In this section we compare the evolution of phase-space DFs for all 
four of the DM halos described in §[2l 

Figure[7]compares the volume distribution of phase space den- 
sity V{f) for the four different halos at z = 0. All three halos from 
the L25 simulation (Gl, G2, G3) show almost identical profiles 
in V{f) confirming the universality of the process that produced 
the phase-space DF. For halo L20, V[f) lies systematically above 



the other curves especially at higher values of /, where it also ex- 
tends to large values of /. This is numerical consequence of the 
increased mass resolution o f the simulation as was first shown by 
ISharma & Steinmet3 j2006l) . This indicates that the absolute value 
of / derived in the previous section is somewhat dependent on the 
mass resolution of the simulation and increases slightly with in- 
creasing mass resolution. In all four halos V(f) is well approx- 
imated by a power-law over nearly 6 orders of magnitude in /. 
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Figure 7. V{f) for four halos, at 2 = 0. The triple-dot-dashed Une is for 
Gl, the dotted line for G2, the dashed line for G3, and the solid line is for 
higher mass resolution simulation L20. 



Figure 8. Q (thin lines) and F (thick Unes) for all the halos studied, at 
z = 0. Triple-dot-dashed lines are for Gl, dotted lines for G2, dashed lines 
for G3, and the solid lines are for L20. 



The first column of Table [T] gives details of the power-law fits to 
V{f) (i.e. the power-law slopes and their errors over the range 
10"* < / < 10^-^.) The power-law sl opes we obtained are similar 
to those obtained by previous authors jArad et aljiooj) . 

In Figure[8]we plot F{r) and Q{r) for the four different halos 
ai z — 0. The top set of four curves show Q{r) while the lower 
set of curves show F{r). In all four halos Q{r), is well fitted by a 
power-law of slope /? = 1.8 — 1.9 (see Table[TJ. F{r) shows signif- 
icant deviations from a simple power-law profile, with a systematic 
upturn beyond r/r^ir > 0.1. The higher resolution simulation L20 
appears to have systematically higher F and Q values than the ha- 
los from the low resolution simulation (a numerical conse quence of 
the higher mass resolution. ISharma & Steinmetj ( l2006h .) Table [T] 
gives the values for the slopes and the error-bars on the power-law 
fits to Q{r) and F{r) for r < 0.6rvir as well as the inner power- 
law slope of _F(r-). 



5 ON THE INTERPRETATION OF Q{R) AND 
PHASE-SPACE DENSITY DISTRIBUTIONS 

Several previous studies have attempted to account for the ori- 
gin of the power-law Q{r) profiles of DM halos. Notably, it has 
been argued that power-law profiles result from virialization and 
not from the hierarchical sequence of mergers, since they are also 

i iroduced in simple spherical gravitationa l collapse simulations 
Tavlor & Navarrol200ll:[Bame s et al. 2006. 2007). Our results pre- 
sented in the previous section (Figure|5j show that while Q{r) and 
F{r) have approximately power-law form within 0.6rvir at a given 



redshift, F{r) flattens out and remains quite flat beyond this ra- 
dius. Furthermore, as the hierarchical growth of the halo progresses 
these approximately power-law profiles extend to larger radii until 
they encompass all the mass within the virial radius at z = 0. In 
this section we attempt to understand the origin of approximately 
power-law profiles. 

It has been sho wn fro m both theoretical arguments and nu- 
merical simulations jPehn en 2005; Vass et al. 2008) that power- 
law profiles of phase-space density for central cusps are well pre- 
served during major mergers. At intermediate times during a major 
merger, deviations from the power-law profile are seen but these 
largely disappear at the end of the merger. In major mergers there 
are two main reasons for the preservation of the power-law profiles. 
First, the most tightly bound material - that forming a steep central 
cusp or shallow core preserves its phase-space density in the final 
remnant. This is a consequence of the additivity of the excess mass 
function - a result of the fact that steeper cusps are less mixed than 
shallower cusps (Dehnen 2005). Over 60% of the material in the 
central cusp (within one scale radius) of a prog enitor NFW halo re - 
mains within the cusp of the merger remnant dValluri et al.ll2007h . 
The second reason for preservation of power-law profiles at large 
radii is that material outside 3 scale radii of the progenitor halos is 
redistributed to other radial bins almost uniformly from each radial 
interval. In addition nearly 40% of the material within the virial ra- 
dius of the progen itor halos is ejected to beyond the virial rad ius of 
the final remnant jKazantzidis et alJl2006l : IValluri et al" ]|2003), and 
this material can be shown to have originated in roughly equal frac- 
tions from each radial interval beyond 3 scale radii and has higher 
phase-space density than expected from the simple power-law ex- 
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Table 1. Power-law indices for fits to V, Q, and F for four halos at z = for r < O.Grvir 



Halo 


V 


Q 


F 


Gi 


-2.34 ±0.02 


-1.84 ±0.01 


-1.59 ± 0.05 


G2 


-2.37 ±0.02 


-1.82 ±0.01 


-1.46 ± 0.06 


Ga 


-2.35 ±0.02 


-1.75 ±0.01 


-1.42 ± 0.03 


L20 


-2.27 ±0.01 


-1.87 ±0.01 


-1.64 ±0.07 



trapolation of the inner power-law. This self-similar redistribution 
of material contributes to the preservation of power-law profiles in 
Q{r) and F{r) in equal mass binary mergers. Thus, once a power- 
law phase-space DF has been established in a DM halo major merg- 
ers will not destroy this profile. 

As was discussed in the previous section (Figure [SJl, the 
coarse-grained phase-space density / in ACDM halos has a log- 
normal distribution with the median of log(/) (the peak of the 
histogram /peak) evolving steadily toward lower values of / as 
the halo grows. This evolution results from both the decrease in 
/ due to mixing when matter accretes onto a halo and virializes 
and due to the decrease in the mean density as the Universe ex- 
pands (cosmic expansion reduces the mean density of the Universe 
as p{z) oc (1 + z^.) Decrease in the median phase-space density 
/peak with decreasing redshift should therefore be at least as fast as 
the decrease in the mean density of the Universe, or faster. The de- 
crease of phase-space density / due to mixing is the consequence 
of complex non-linear evolution, hence it is difficult to arrive at a 
analytic expression to represent the cosmological evolution of the 
median phase-space density /peak- We therefore attempt to derive 
it empirically. 

Figure |9] shows the median phase-space density /peak at each 
redshift derived from the histograms of log(/) in Figure|6]for each 
of the 4 Milky-Way sized cosmological halos in this study. The 
points represent the values of /peak as a function of expansion fac- 
tor a=(l+2)~^, while the thin dot-dash lines connect points for 
each of the 4 halos. The thick solid line is the best fit power-law for 
all halos taken together: /peak(ffl) oc a^* "^*^'^. As expected, this 
decrease is faster than a~'^. 

Thus on average, as a halo grows via accretion its median 
density decreases as a power-law with time, despite the fact that 
the most centrally concentrated material retains its original high 
phase-space density. This power-law profile in /peak leads to some 
insights on the development of phase-space density pro files. We 
can b reak down the formation of halos into two phases ihi et al.l 
l2007h : the fast accretion regime during which halo mass grows very 
rapidly and the slow accretion regime. In the fast accretion regime, 
the mixing processes are very efficient as the potential well is estab- 
lished and potential fluctuates rapidly and constantly. This results 
in a rapid decrease in the over all central phase-space density of 
the halo (as seen from the rapid drop in the central most regions of 
the profiles in Fig.|2l) The inner profile of phase-space distribution 
should be largely set during this stage. As halos grow subsequently. 




Figure 9. The median phase-space density (peak of histograms in Fig. |6) 
(/peak) 3t each redshift as a function of the expansion parameter a = (1 + 
z)^^ for four Milky-way sized halos (points). The thin dot-dashed lines 
connect points for individual halos. The best-fit (solid line) has a slope of 
-4.3 ± 1.1. 

either by major mergers or quiescent accretion of smaller halos, the 
average phase-space density decreases, but the pre-existing high 
central phase-space density cusps are are preserved. 

At these later times the evolution of the true phase-space den- 
sity is complex and occurs due to the accretion of high phase-space 
density material in subhalos as well losely bound material at the 
edges of the halo. The steady decrease in the amplitude and increase 
in the slope of the F{r) profiles in Fig. |2] with redshift show that 
the true phase-space density does not obey the power-law profiles 
seen in Q(r) at all redshifts. Variations in the individual power-law 
profiles of halos both at z = (Fig. [8} and with redshift reflect the 
large variation in the cosmic accretion histories of individual halos. 
Additional deviation is due to the matter bound to subhalos that 
survives with high phase-space density and leads to the variation of 
the coarse-grained / of some six orders of magnitude in the outer 
regions of halos. 

All this indicates that the actual phase-space distribution is 
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not as universal and simple as the Q{r) profiles lead one to believe. 
The mechanism behind the universal power-law form of the Q{r) 
profiles is therefore likely to be different (e.g., it is not affected by 
the presence of subhalos) and simpler than the processes that shape 
the distribution of the coarse-grained phase-space density. The Q 
profile is a ratio of the density and velocity dispersion (which can 
be interpreted as a measure of temperature) and is therefore related 
to the entropy. For an ideal monatomic gas the entropy can be de- 
fined as ifgas = T/ f?^'^ , while for dark matter, by analogy, the en- 
tropy can be defin ed as j't'dm = cP ' j j?J^ jpaltenbacher et al.l2007h . 
This, as noted bvl Hoffman etalj (12007'), gives T^dm « Q"^/^ or 
Katci oc r^'^ for Q oc r~^'*. This power-law form and slope of 
the entropy profile is very similajf] to the one found for the gas 
in the outer regions of clusters in co smological simulations (e.g., 
iBorgani et al]|2004 IVoit et al.ll20o3) . This power-law is also pre- 
dicted by the models of spherical accretion jTozzi & Normanl200lh 
and reflects the encreasing entropy to which the accreting material 
is heated as the halo grows its mass. Although the processes gov- 
erning the virialization of gas and dark matter are different (the 
short-range local ineractions for the former, and long-scale interac- 
tions for the latter), the fact that the resulting entropy profiles are 
quite similar indicate that they lead to the same distribution of en- 
tropy. The Q(r) profile therefore may reflects the overall entropy 
profile of dark matter, not the coarse-grained local phase-space den- 
sity, which exhibits a more complicated behavior. 



6 SUMMARY AND CONCLUSIONS 

We have investigated the evolution of the phase-space density of 
the dark matter in cosmological simulations of the formation of 
Milky- Way- sized dark matter halos. The analysis was carried out 
using two different codes for estimating the phase-space density. 
Both codes give qualitatively similar results, but the estimated val- 
ues of phase- space density / are quite sensitive to the type of code 
and, for a given code, also depend quite sensitively on the choice of 
smoothing kernel used. Based on comparisons of the two codes (Fi- 
EstAS and EnBiD ) and various smoothing parameters (Appendix), 
we select the EnBiD code with n = IQ kernel smoothing and 
present results for the analysis with this set of parameters. 

The simulations presented in this paper complement our anal- 
ysis of the evolution of phase-space density in binary major merg- 
ers jVass et al] '2008). We confirm that the profiles of Q{r) = 
Pdm/o"dm computed by previous authors can be described by a 
power-law Q{r) oc i over more than two orders of mag- 

nitude in radius in all halos. The median of the phase-space den- 
sity (F{r)) at given radius r, however, exhibits a more compli- 
cated behavior. Although F{r) is approximately a power-law for 

^ The slope is even more similar if on e takes into account the ra ndom bulk 
motions of the gas in estimating Kg^s jpaltenbacher et alj2007l) . 



r < 0.6rvir, the profiles generally flatten in the outer regions. 
Subhalos contribute somewhat to this behavior, although their ef- 
fect is limited by the relatively small fraction of mass (< 0.1) 
bound to them. However, in addition to subhalos, a significant 
fraction of high phase-space density matter is in the relatively un- 
mixed streams of alre ady disrupted subhalos (e.g.. lArad et al.l2004l ; 
lOiemand et alj|20o3) . The fraction of mass in such streams can be 
substantial at large radii. 

This behavior holds at earlier epochs. From z — 5 to z = 0, 
material within rvii{z) at each redshift follows a power-law in Q 
with an approximate power-law slope of ~ —1.8 to —1.9. In con- 
trast, F{r) can only be well described by a power-law in the in- 
ner regions, and its slope changes continuously with redshift. Be- 
yond the virial radius, Q (a quantity that is obtained by averaging 
over increasingly large volumes) decreases rapidly with radius, but 
the median value of F flattens significantly. We argue that F is a 
more physically meaningful quantity, especially for understanding 
the evolution of phase-space density in collisionless DM halos, as 
it measures the median of the true coarse-grained phase-space den- 
sity. 

At all redshifts, the highest values of phase-space density / 
are found at the centers of dark matter subhalos. In the center 
of the main halo, the median phase-space density (F) drops by 
about an order of magnitude from z = 9 to 2 = 0. In con- 
trast, the centers of DM subhalos maintain their high values of 
/ ~ lO^M0kpc"^(kms"^)"^ at aU redshifts. The highest val- 
ues of / at the center of the main halo are, therefore, lower and less 
representative of the primordial phase-space density of dark matter 
particles than the central value of / in the high phase-space den- 
sity subhalos. At r^ir, the decrease in median phase-space density 
is much more significant, with F ~ 3OM0kpc~^(kms~^)~^ at 
z = 9 decreasing to _F « lO~^M0kpc~^(kms"^)"^ at z = 0, a 
decrease of over four orders of magnitude (Figure|4ll. 

The evolution of F{r) and V{f) with redshift are consistent 
with expectations from the Mixing Theorems, which require that 
mixing reduces the overall phase-space density of matter in col- 
lisionless systems and that the volume of phase-space associated 
with any value of / increases due to mixing and relaxation. 

The distribution of / is approximately log-normal until 2 ~ 3. 
As time progresses, the mean and median of log(/) shift to pro- 
gressively lower values as a larger and larger fraction of matter 
undergoes mixing and moves to lower values of /. Some frac- 
tion of high phase-space density material does survive in the cen- 
ters of subhalos and in the relatively unmixed streams leftover 
after subhalo disruptions, which skews the distribution. Remark- 
ably, the highest phase-space density particles at z — have re- 
tained their phase-space density since z ~ 9, the earliest epoch 
we analyzed. The phase-space density in the centers of dark mat- 
ter subhalos is therefore representative of the mean phase-space 
density of DM at the highest redshifts. This can potentially al- 
low for stronger constraints to be placed on the nature of DM 
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particles from the Tremain -Gunn bound jTremaine & Gunnll 19791 : 
ikogan & DalcantorJbOOOl) . 

The median value of phase-space density /peak decreases with 
decreasing redshift and is given by the power-law /peak (a) oc 
^-4.3±i.i gjj^ jj^g result of two different processes; (a) mixing 
within virialized halos which reduces the coarse-grained phase- 
space density of matter that has turned around from the Hubble 
flow and (b) the cosmic expansion of the Universe which results in 
a decrease in the configuration-space density of "unprocessed" ma- 
terial (material that has yet to fall into a massive halo and become 
virialized). Thus, as a halo grows by accretion and the Universe ex- 
pands, the median phase-space density in halos decreases steadily. 
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APPENDIX A: COMPARISON OF "FiEstAS " AND 
"EnBiD " ANALYSIS OF HALO Gl 

The numerical estimation of coarse-grained phase-space densities 
during the evolution of the 4 N-body halos presented in this pa- 

S ier was carried out using tw o publicly availa ble codes "FiEstAS " 
Ascasibar & BinnevI l2005h and "EnBiD " l lSharma & Steinmetj 
1200^ . A comparison of these two codes has previously been pre- 
sented by Shcurma & Steinmetz (2006) who showed that for an an- 
alytic DF (e.g. for the spherical Hemquist potential), that "En- 
BiD with n = 10 kernel smoothing " gave the highes t fidelity to 
the analytical DF. In a related work Vass et all j2008l) confirmed 
the findings of Sharma & Steinmetj ^2006) for a spherical isotropic 
NEW halo. This latter study is the main basis for our choice of "En- 
BiD with n = 10 kernel smoothing" . 



Table Al. Power-law indices V{f) oc / " and F(r) oc r ^ at z 
from different codes 



Code 


a 


/3 


FiEstAS 


2.62 ±0.06 


1.43 ± 0.02 


EnBiD (no smoothing) 


2.35 ±0.02 


1.65 ± 0.02 


EnBiD (FiEstAS smoothing) 


2.46 ± 0.03 


1.69 ± 0.07 


EnBiD (kernel n = 10) 


2.34 ±0.02 


1.59 ± 0.05 



However, it is unclear whether the comparisons with analytic 
profiles of isolated halos at z — are valid for matter distributions 
arising from cosmological N-body simulations, at high redshifts, 
where the majority of the particles actually lie outside virialized 
halos. Our purpose in this Appendix is to present a comparison 
of results obtained with the different codes at a range of redshifts 
to allow readers to appreciate how sensitive some of the results 
presented in this paper are to the choice of code and smoothing 
parameters used in density estimation. 

Since the coarse-grained DF /(x, v) is a 6-dimensional func- 
tion, it is difficult to compare estimates for this function obtained 
using different codes. Traditionally the best single variable func- 
tion to compare is the volume DF V{f). Variations in the estima- 
tion of / translate to variations in estimation ofV{f). In Figure IaTI 
we compare the volume density of phase-space V{f) obtained for 
the Gl halo at 4 different redshifts using the FiEstAS code (triple- 
dot-dashed curves), EnBiDcode with no smoothing (dot-dashed 
curves), EnBiD with FiEstAS smoothing (long-dashed curves), and 
EnBiD with a n = 10 kernel ( dotted curves). In each case 
we see that V{f) has a nearly power-law distribution over more 
than six orders of magnitude in / at 2: — (from ~ 10~*- 
lO^A/0kpc"^(kms"^)"^). Table IaT] gives the slopes of power- 
law fits for the four different estimates of V{f) over the range 
10"'' < / < 10^ at z = 0. The FiEstAS estimate of V{f) is sys- 
tematically lower than all EnBiD estimates at high values of / (at 
all redshifts). In addition all EnBiD curves extend to much higher 
values of / than the FiEstAS estimate (this is particularly true at 
z = 9 where the FiEstAS estimate differs from the other estimates 
both quantitatively and qualitatively). The results from the vari- 
ous EnBiD estimates differ very little at intermediate values of / 
(10~^-10^) at 2 = and consequently we are confident that con- 
clusions drawn from the median / is quite insensitive to the details 
of the EnBiD parameters used to obtain /. 

The differences between the various estimates of F at large 
radii become signficantly larger at higher redshift. In Figure lA2l we 
plot the four different estimates of F(r) at 4 different redshifts in 
the evolution. The vertical dashed line in each panel represents the 
virial radius of the main halo at that redshift. We see that F from 
any of the codes is quite flat beyond the virial radius in all cases, 
but the absolute values of the curves differ significantly. Note that 
at 2 = 9, the different estimates can differ by as much as two orders 
of magnitude. 
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Figure Al. The volume DF y (/), at different redshifts during the evolution 
of the Gl halo in the L25 run, resulting from the FiEstAS and EnBiD with 
three different parameters as indicated in the line legends. 
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Figure A2. F for the Gl halo in the L25 simulation at different redshifts 
obtained using four different analyses as indicated by line-legends. The 
dashed vertical line is situated at the vmal radius of the virialized main 
halo at each redshift. 



In this paper we choose to present results from EnBiD with 
n = 10 kernel smoothing over the other estim ates largely because it 
does the best job of reproducing analytic DFs fe harma & Steinmetd 
I2OO6) and because it appears to provide a good upper limit to the 
phase-space density both at high and low values of /. In the ab- 
sence of an analytic comparison of the estimates at high redshift, 
we caution the reader to refrain from drawing very strong conclu- 
sions regarding absolute values of / or F from the results presented 
here. 
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